AirTemperatureDailyMeanRead Subroutine

SUBROUTINE AirTemperatureDailyMeanRead &
!
( time, dem )

IMPLICIT NONE

!Arguments with intent(in):
TYPE (DateTime), INTENT(IN) :: time !!current time
TYPE (grid_real), INTENT(IN) :: dem  !digital elevation model to apply drift

!local declarations:
TYPE (DateTime) :: time_toread !! time to start reading from
CHARACTER (LEN = 300) :: filename
CHARACTER (LEN = 300) :: varname
REAL (KIND = float)   :: rsquare
INTEGER (KIND = short) :: i, j

!-------------------------end of declarations----------------------------------

IF ( .NOT. (time < timeNew) ) THEN
   
   !time_toread = time + - (dtTemperatureDailyMean - thermometersDailyMean % timeIncrement)
   time_toread = time
   timeString = time_toread
   CALL Catch ('info', 'AirTemperatureDailyMean',   &
		                 'read new temperature data: ', &
                     argument = timeString)
   
   !update lapse rate from maps when required
    IF (elevationDrift == 1) THEN
       IF (time == lapse_map % next_time) THEN !update lapse_map
          varname = lapse_map %var_name
          filename = lapse_map % file_name
          !destroy current grid
          CALL GridDestroy (lapse_map)
          !read new grid in netcdf file
          CALL NewGrid (layer = lapse_map, fileName = TRIM(filename), &
                         fileFormat = NET_CDF, &
                         variable = TRIM(varname), &
                         time = time_toread)
      END IF
    END IF

   SELECT CASE (interpolationMethod)
    CASE (0) !input grid
      
        CALL ReadField (fileGrid,  time_toread, &
                        dtTemperatureDailyMean, dtGrid, &
                        'M', temperatureAirDailyMean, &
                         varName = temperatureAirDailyMean % var_name)
 
    CASE DEFAULT !requires interpolation
      !read new  data
      CALL ReadData (network = thermometersDailyMean, fileunit = fileunit, &
                      time = time_toread, aggr_time = dtTemperatureDailyMean, &
                      aggr_type = 'ave', tresh = valid_prcn)
      
      IF (elevationDrift == 1) THEN
        !update lapse rate from data  
        IF (lapse_rate_computation == 1) THEN
            
          DO i = 1, thermometersDailyMean % countObs
               vectorT (i) = thermometersDailyMean % obs (i) % value
               vectorZ (i) = thermometersDailyMean % obs (i) % xyz % elevation
          END DO
          
          lapse_computed = LinearRegressionSlope (x=vectorZ, y=vectorT, &
                    nodata = thermometersDailyMean % nodata, r2 = rsquare)
          IF (rsquare < r2min) THEN
              lapse_computed = lapse
          ENDIF
          
          lapse_map = lapse_computed
        END IF

        !shift  observation to reference elevation by applying lapse rate
        CALL ShiftMeteoWithLapse (thermometersDailyMean, lapse_map, refElevation, &
                                  stationsRefElev)
        
        !interpolate 
        IF (interpolationMethod_assignment == 1) THEN !only one method is applied
                
                CALL Interpolate (network = stationsRefElev, &
                            method = interpolationMethod, &
                            near = neighbors, &
                            idw_power = idw_power, & 
                            anisotropy = krige_anisotropy, &
                            varmodel = krige_varmodel, &
                            lags = krige_lags, &
                            maxlag = krige_maxlag, &
                            grid = temperatureAirDailyMean, &
                            devst = grid_devst)                  
                
        ELSE 
            !loop trough interpolation methods
            DO i = 1, 3
                IF (interpolationMethod_vector(i) > 0) THEN
                      
                    CALL Interpolate (network = stationsRefElev, &
                            method = interpolationMethod_vector(i), &
                            near = neighbors, &
                            idw_power = idw_power, & 
                            anisotropy = krige_anisotropy, &
                            varmodel = krige_varmodel, &
                            lags = krige_lags, &
                            maxlag = krige_maxlag, &
                            grid = interpolatedMap (i), &
                            devst = grid_devst) 
                      
                END IF
            END DO

            !compose the final interpolated field
            DO i = 1, interpolationMethod_map % idim
                DO j = 1, interpolationMethod_map % jdim
                    IF (interpolationMethod_map % mat(i,j) /= &
                        interpolationMethod_map % nodata ) THEN
                        temperatureAirDailyMean % mat (i,j) = &
                        interpolatedMap (interpolationMethod_map % mat(i,j)) % mat(i,j)
                    END IF
                END DO
            END DO
        END IF

        !Shift back interpolated field to terrain surface
        CALL ShiftBackWithLapse (temperatureAirDailyMean, dem, &
                                lapse_map, refElevation)
       
         
      ELSE !elevationdrift = 0
        

        IF (interpolationMethod_assignment == 1) THEN !only one method is applied
                
                CALL Interpolate (network = thermometersDailyMean, &
                            method = interpolationMethod, &
                            near = neighbors, &
                            idw_power = idw_power, & 
                            anisotropy = krige_anisotropy, &
                            varmodel = krige_varmodel, &
                            lags = krige_lags, &
                            maxlag = krige_maxlag, &
                            grid = temperatureAirDailyMean, &
                            devst = grid_devst) 
                
        ELSE
            !loop trough interpolation methods
            DO i = 1, 3
                IF (interpolationMethod_vector(i) > 0) THEN
                      
                    CALL Interpolate (network = thermometersDailyMean, &
                            method = interpolationMethod_vector(i), &
                            near = neighbors, &
                            idw_power = idw_power, & 
                            anisotropy = krige_anisotropy, &
                            varmodel = krige_varmodel, &
                            lags = krige_lags, &
                            maxlag = krige_maxlag, &
                            grid =  interpolatedMap (i), &
                            devst = grid_devst) 
                    
                END IF
            END DO

            !compose the final interpolated field
            DO i = 1, interpolationMethod_map % idim
                DO j = 1, interpolationMethod_map % jdim
                    IF (interpolationMethod_map % mat(i,j) /= &
                        interpolationMethod_map % nodata ) THEN
                        temperatureAirDailyMean % mat (i,j) = &
                        interpolatedMap (interpolationMethod_map % mat(i,j)) % mat(i,j)
                    END IF
                END DO
            END DO

        END IF  !1 or more interpolation methods       
      END IF  !elevationDrift
		
	END SELECT

  !apply scale factor and offset, if defined
	IF (offset_value /= MISSING_DEF_REAL) THEN
	   CALL Offset (temperatureAirDailyMean, offset_value)
	END IF
	
	
	IF (scale_factor /= MISSING_DEF_REAL) THEN
	   CALL Scale (temperatureAirDailyMean, scale_factor)
	END IF


  !grid exporting
  IF(export > 0 ) THEN
	  IF( time == timeNewExport .AND. &
        time >= export_start .AND. &
        time <= export_stop ) THEN
        timeString = time
        timeString = timeString(1:19)
        timeString(14:14) = '-'
		    timeString(17:17) = '-'
        
        grid_devst % reference_time = temperatureAirDailyMean % reference_time
        IF (needConversion) THEN
           CALL GridConvert (temperatureAirDailyMean, exportedGrid)
           CALL GridConvert (grid_devst, exportedGridVar)
        ELSE
           exportedGrid = temperatureAirDailyMean
           exportedGridVar = grid_devst
        END IF 

        SELECT CASE (export_format)
        CASE (1) !esri-ascii
              export_file = TRIM(export_path) //  TRIM(timeString) // &
                           '_temperature_daily_mean.asc'
              CALL Catch ('info', 'AirTemperatureDailyMean',   &
		                       'exporting grid to file: ', &
                           argument = TRIM(export_file))
              CALL ExportGrid (exportedGrid, export_file, ESRI_ASCII)
              
              IF (krige_var == 1) THEN !export kriging variance
                    export_file_var = TRIM(export_path) //  TRIM(timeString) // &
                           '_temperature_daily_mean_variance.asc'
                    CALL Catch ('info', 'AirTemperatureDailyMean',   &
		                       'exporting variance grid to file: ', &
                              argument = TRIM(export_file_var))
                    CALL ExportGrid (exportedGridVar, export_file_var, ESRI_ASCII)
              END IF
              
        CASE (2) !esri-binary
              export_file = TRIM(export_path) //  TRIM(timeString) // &
                           '_temperature_daily_mean'
              CALL Catch ('info', 'AirTemperatureDailyMean',   &
		                       'exporting grid to file: ', &
                           argument = TRIM(export_file))
              CALL ExportGrid (exportedGrid, export_file, ESRI_BINARY)
              
              IF (krige_var == 1) THEN !export kriging variance
                   export_file_var = TRIM(export_path) //  TRIM(timeString) // &
                               '_temperature_daily_mean_variance'
                   CALL Catch ('info', 'AirTemperatureDailyMean',   &
		                       'exporting variance grid to file: ', &
                              argument = TRIM(export_file_var))
                   CALL ExportGrid (exportedGridVar, export_file_var, ESRI_BINARY)
              END IF
              
        CASE (3) !net_cdf
              CALL SetCurrentTime (time, exportedGrid)
              CALL Catch ('info', 'AirTemperatureDailyMean',   &
                           'exporting grid to file: ', &
                           argument = TRIM(export_file))
              CALL ExportGrid (exportedGrid, export_file, 'air_temperature', 'append')
              
              IF (krige_var == 1) THEN !export kriging variance
                  CALL SetCurrentTime (time, exportedGridVar)
                  CALL Catch ('info', 'AirTemperatureDailyMean',   &
		                       'exporting grid to file: ', &
                              argument = TRIM(export_file_var))
                 CALL ExportGrid (exportedGridVar, export_file_var, 'air_temperature_variance', 'append')
              END IF
        END SELECT
        timeNewExport = timeNewExport + export_dt
    END IF
  ENDIF



  !time forward
  timeNew = timeNew + dtTemperatureDailyMean
END IF

RETURN
END SUBROUTINE AirTemperatureDailyMeanRead